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Abstract 

Many boundary element integral equation kernels are based on the Green’s functions of the Laplace and Helmholtz 
equations in three dimensions. These include, for example, the Laplace, Helmholtz, elasticity, Stokes, and Maxwell’s 
equations. Integral equation formulations lead to more compact, but dense linear systems. These dense systems are 
often solved iteratively via Krylov subspace methods, which may be accelerated via the fast multipole method. There 
are advantages to Galerkin formulations for such integral equations, as they treat problems associated with kernel 
singularity, and lead to symmetric and better conditioned matrices. However, the Galerkin method requires each entry 
in the system matrix to be created via the computation of a double surface integral over one or more pairs of triangles. 
There are a number of semi-analytical methods to treat these integrals, which all have some issues, and are discussed 
in this paper. We present novel methods to compute all the integrals that arise in Galerkin formulations involving 
kernels based on the Laplace and Helmholtz Green’s functions to any specified accuracy. Integrals involving completely 
geometrically separated triangles are non-singular and are computed using a technique based on spherical harmonics 
and multipole expansions and translations, which results in the integration of polynomial functions over the triangles. 
Integrals involving cases where the triangles have common vertices, edges, or are coincident are treated via scaling and 
symmetry arguments, combined with automatic recursive geometric decomposition of the integrals. Example results 
are presented, and the developed software is available as open source. 


1 Introduction 

The Galerkin boundary element method (BEM) is a powerful method for solving integral equations involving kernels 
based on the Laplace and/or Helmholtz equations’ Green’s functions in three dimensions nil Em. When the boundary 
is discretized using triangular elements, constructing the system matrix requires computing double surface integrals over 
pairs of these triangles. Because the kernels being integrated are singular, these integrals can be difficult to compute, 
especially when the two triangles are proximate, share a vertex, an edge, or are the same. Depending on the relative 
geometry of the two triangles, there are many different methods for computing them. Eor example, when the two 
triangles do not touch, the integral is completely regular and can be computed accurately via numerical means, e.g., 
Gaussian quadrature H. Semi-analytical methods, where the inside integral is computed analytically and the outside 
integral is computed numerically, have also been proposed lEia. 

However, when the two triangles share a vertex, an edge, or are the same, things become more complicated. There 
are analytical expressions for the case when the two triangles are the same (71, but not for when they share only a vertex 
or an edge. In these cases, the semi-analytical methods do not always work. This is because, depending on the kernel 
being integrated, the inside integral can be hypersingular. While there are analytical expressions available for them, 
they are singular along the corners and edges of the corresponding triangle. When the two triangles share a vertex or an 
edge, these singularities are included in the outside integral. The usual semi-analytical methods will not work in these 
cases because they are not designed to properly handle the singularities. 
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Figure 1. Several different types of boundaries: Fi is closed; r 2 is closed, but contains another boundary. Fa, also 
closed; F 4 and F 5 form a closed region, but this region is divided into two by Fg; and Fy is open. 


The double integrals are weakly singular, so while the inside integrals may be hypersingular and the expressions 
for them may be singular in some places, they are completely integrable. Nevertheless, actually integrating them in 
practice can be hard. Therefore, more sophisticated semi-analytical methods have been developed over the years. These 
include: singularity subtraction and “to the boundary” techniques lISlI^ fTOlfTTlfT^ : specialized quadrature methods 
that are designed for the singularities involved, such as those based on the double exponential formula imiia; and 
other regularization methods, such as the Duffy transformation ifTSlflblfTTll . Many of these methods work very well, 
but because they all attempt to tackle the singularity issue directly, their analysis is very involved. 

In this paper, we present a method for computing the integrals that completely avoids the computation of singular 
integrals. The approach relies on several scaling properties of the integrals and the kernels being integrated. When two 
triangles share a vertex, an edge, or are the same, the integral is decomposed into several smaller integrals, some of 
which are related back to the original integral using simple analysis. This is done in such a way so that only regular 
integrals need to be computed explicitly. Any integrals involving singularities are computed implicitly during the 
procedure. The regular integrals can be computed using standard semi-analytical methods, but in this paper, we also 
present an analytical method for doing so. This method uses spherical harmonics and multipole and local expansions 
and translations. The only source of error in the method is from truncating these expansions. However, this error is 
precisely controlled by choosing the appropriate truncation number or recursively subdiving the problem. Although we 
developed these methods primarily for kernels related to the Laplace equation’s Green’s function, we also describe how 
they can be extended to kernels related to the Helmholtz equation’s Green’s function. Finally, we have implemented 
a Galerkin BEM library in MATLAB for the Laplace equation using these two methods. We show some example 
problems and provide some error analysis. We have made this library freely available for download ifTSl . The library 
can be used to recreate all the examples seen in this paper or to create entirely new ones. 

2 Background 

The Galerkin BEM is a powerful method for solving integral equations related to the Laplace or Helmholtz Green’s 
functions in three dimensions. Consider the following boundary value problem (BVP): 

V^(/>(x) = 0 or V^(/> (x)-I-(x) = 0, X e H, (1) 

where A: > 0 is the wavenumber in the Helmholtz equation. General boundary conditions are enforced on both sides of 
the boundary, and are given by: 

Oil (x) (x) -f (x) 9 + (x) -f aj" (x) (j}~ (x) -f (x) q~ (x) = 71 (x) , x G F, (2) 
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0^2 (x) (x) + /3^ (x) 9 + (x) + a 2 (x) (j) (x) + (x) q (x) = 72 (x) , x S F. (3) 

In addition, boundary conditions must be set at infinity. For the Laplace equation, the potential should decay to zero at 
large distances. For the Helmholtz equation, the potential should decay to zero at large distances as well, but should 
also be composed of outgoing waves only. The Sommerfeld radiation condition provides such a constraint; 


lim (p (x) = 0 , 

|x|—>-oo 


lim |x| 

|x|—>-oo 



ik(f> (x) I = 0. 


(4) 


The boundaries can be closed or open. See Fig.[T]for some example boundaries. For closed boundaries, the interior 
side of the boundary is the ” side and the exterior side is the “+” side. For open boundaries, since there is no inside 
or outside, designating each side of the boundary as “+” or ” can be done arbitrarily. The values, and (j)~, are 
the potential on the “+” and ” sides of the surface, respectively. Likewise, the values, 9 + and q~, are the normal 
derivatives of the potential on the “+” and ” sides of the surface, respectively, and are given by 


9 


± 


d(j)^ 

dn^ 


(n±.Vx) 


(5) 


where = =Fn (i.e., n''' goes from the “+” side to the ” side, and vice versa). 

To solve the BVP, we use an indirect boundary integral formulation called the layer potential formuation. Using 
Green’s theorem, the Laplace or Helmholtz equation is transformed from a differential equation into an integral equation: 


0 (x) = L [cr] (x) + M [^J\ (x), 


( 6 ) 


where 


L [cr] (x) = / cr (x') G {x — x') dS (x'), 
Jx'GF 

M [/r] (x) = f /c (x') (n'• Vx') G (x — x') els'(x') 
Jx'er 


are the single- and double-layer potentials 11911 . and 


C(r) = ^ or c(r) = '■■I* 


47r r 


Att r 


(7) 

( 8 ) 

(9) 


is the Laplace or Helmholtz equation’s Green’s function. The single-layer potential, L [a] (x), is the potential due to the 
monopole source density distribution, cr (x'), on the boundary. Likewise, the double-layer potential, M [/i] (x), is the 
potential due to the dipole source density distribution, /i (x'), on the boundary. In the differential equation, we seek 
a solution to the potential governed by the Laplace equation. However, in the integral equation, we seek the source 
density distributions, cr (x') and /c (x'), on the boundary that give rise to that potential. The advantage of the BEM is 
that the expression in Eq. (|^ relating the source density distributions back to the potential automatically satisfies the 
original differential equation. Moreover, the Green’s functions satisfy the boundary conditions at infinity, so as long as 
the source density distributions are bounded and finite, the single- and double-layer potentials will satsify them as well. 
Thus, we need only concern ourselves with searching for the source density distributions that satisfy the remaining 
boundary conditions. To do this, we need to express the potentials and normal derivatives on either side of the boundary 
in terms of the source density distributions. Jump conditions provide such a relationship; 

(j)^ (x) = L [cr] (x) -b M [n] (x) ± (x), (10) 

q^ (x) = tL' [cr] (x) =F M' [/r] (x) -b (x), (11) 

where 

L'[cr] (x) = (n • Vx) / (j{x')G{x-x')dS{x'), (12) 

Jx'GF 

M'[/r] (x) = (n • Vx) / /r (x') (n'• Vx') G (x - x') dS” (x'). (13) 

dx'Gr 
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Plugging Eqs. ( [T0| and (111 into the boundary conditions and rearranging, 

fli (L [a] + M [^]) + bi {L' [a] + M' [/r]) + citr + difi = 71 , (14) 

02 {L [a] + M [^]) + 62 (U' [a] + M' [/r]) + C 2 cr + d 2 H = 72, (15) 

where 

ai=a^+a^, 6 i = -/3++/3f, ci = -{p^+P^), di = -{at-a^), (16) 

a 2 = a^+a^, 62 = -/3^+/3^, C2 = ^ + P 2 ) ^ ^(2 = ^ • (17) 

In these expressions, we have dropped the argument, x, to save space (i.e., a (x) becomes a). In order to make the 
problem computationally tractable, the source density distributions, a (x') and fj. (x'), are each written as a linear 
combination of N basis functions; 


N N 

(^') = ^ (x') • 

i=i i=i 


(18) 


For constant triangular elements, there is one basis function per triangle that is equal to one on that element and zero 
everywhere else. For linear triangular elements, there is one basis function per vertex that is equal to one at that vertex, 
zero at all the other vertices, and piecewise linear everywhere else. We need to compute the coefficients of these basis 
functions so that the boundary conditions are satisfied. In other words, we seek cji, 0 - 2 , ■ ■ •, Qiat, ^ 1 , ^ 2 , • ■ ■, and /tat 
such that 


N 


N 


aj {aiL [fj] + biL' [fj] + afj) + Y N I/f] + [fi\ + c^i) = 7i, 


i=i 

N 


N 


'Y, ^3 [fj] + b2L' [fj] + C2fj) + ^ /Tj {a2M [fj] + b2M' [fj] + ((2) — 72- 
i=i 


Or, more compactly. 


N 


N 


Y = 71 > 


(19) 


( 20 ) 


( 21 ) 


where 


N N 

Y [fj] + Y [fj] = 72 , 

i=i i=i 

^1 [fj] — [fj] + [fj] + 7?! [/_,] = aiM [fj] + biM' [fj] + di, 

^2 [fj] = 0 . 2 L [fj] + b2L' [fj] + C2fj, B 2 [fj] = a2M [fj] + 62 M' [fj] + ^ 2 - 


( 22 ) 

(23) 

(24) 


The two most commonly used methods for enforcing the boundary conditions are the collocation method and the 
Galerkin method. The collocation method works by enforcing the boundary conditions at N matching points: 


N 


N 


Y [fi] [fj] = 71 (x*), i = l,2,...,N, 


j=i 

N 


N 


Y [fj] [fj] (xO = 72 (x,), i = l,2,...,N. 

j=i i=i 


(25) 


(26) 


The method is so named because these points are typically collocated with the modeling elements (e.g., in the case 
of constant triangular elements, one is placed at every triangle’s centroid). Collocation methods have long been used 
(201 EH. They are easy to understand, and the integral expressions necessary for implementing them have been studied 
and derived by many authors. This includes, for example, piecewise constant and linear basis functions on triangular 
elements, which are the most commonly used However, they suffer from a few problems. Many of the 
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(3) (4) 

Figure 2. In practice, the relative geometry of two triangles is one of the following: (1) the two triangles do not touch; 
(2) they share a vertex; (3) they share an edge; or (4) they are the same. 


boundary integrals are hypersingular, which make them hard (or sometimes even impossible) to compute, especially for 
points on the corners or edges of the boundary. 

The Galerkin method overcomes these problems by enforcing the boundary conditions in an integral sense. The 
boundary integral equation is multiplied by each of the same N basis functions from before and integrated over the 
boundary a second time: 


N 

CTJ / fiAi [fj] dS (x) 


N 

i=i 


N 

E [fj] dS (x) 


E^W f^B2[fj]dSi^) 

i=i 



fi-fidS (x), 


fa 2 dS (x), 




(27) 


(28) 


By doing so, all the hypersingular integrals become weakly singular. Morever, the system matrices in the Galerkin 
method are typically symmetric, better conditioned, and have better convergence properties Il26ll27l . However, the 
extra integral over the boundary complicates the computation of the entries in the system matrix. 


3 Double Surface Integrals 

When the boundary is discretized using linear triangular elements, the double surface integrals are performed over pairs 
of these triangles. In each pair, one is called the “source” triangle, and the other the “receiver” triangle. The inside 
integral is over the source triangle, S, and the outside integral is over the receiver triangle, R. Thus, when populating 
the system matrix, we need to compute integrals of the following form: 

I = [ (cto + P • x) / (o-Q + p' • x') F (x - x') dS (x') dS (x), (29) 

Jxgr Jx'gs 

where CTq + p' • x' is the source density distribution over the source triangle, ctq + p • x is the weight function over the 
receiver triangle, and F (r) is the kernel being integrated. To implement the Galerkin BEM described in Sec.|^ we 


need to compute this integral for the following four kernels: 

Fi (x — x') = G (x — x'), (30) 

F 2 (x - x') = (n'• Vx') G(x - x'), (31) 

F 3 (x - x') = (n • Vx) G (x - x'), (32) 

F 4 (x - x') = (n • Vx) (n' • Vx') G (x - x'), (33) 
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Figure 3. For the double surface integral of the Helmholtz equation’s Green’s function, assume the two triangles are 
small and close to each other. In other words, for any two points, one on the source triangle and one on the receiver 
triangle, assume the distance between them, r, is less than 1/fc. 


where Fi (r) and F 2 (r) correspond to the single- and double-layer potentials, and F3 (r) and F 4 (r) correspond to their 
normal derivatives. 

Computing the integral for these four kernels for all commonly encountered geometries is the focus of this paper. In 
practice, the relative geometry of the two triangles is one of the following; (1) the two triangles do not touch; (2) they 
share a vertex; (3) they share an edge; or (4) they are the same (see Fig.[^. These are called the zero-, one-, two-, and 
three-touch cases, respectively. In this naming scheme, the number represents how many vertices the two triangles 
share. For example, in the one-touch case, the two triangles share a single vertex. 


4 Reduction of the Helmholtz Kernel to the Laplace Kernel 


In order to apply the Galerkin BEM to integral equations whose kernels are based on the Helmholtz Green’s function, 
we need to compute double surface integrals of the four kernels given in Eqs. ( |30| l - P3| l over pairs of triangles. The 
remainder of this paper describes several novel methods for computing these integrals when the Green’s function is that 
of the Laplace equation. These methods are not directly applicable to the Helmholtz Green’s function. However, they 
can be extended to the Helmholtz equation (and any other problem involving the Laplace and/or Helmholtz kernels, 
such as elastostatics, Stokes problems, and Maxwell’s equations) via the singularity subtraction technique. 

To reduce the Helmholtz kernel to the Laplace kernel, we use the singularity subtraction technique Eol. As an 
example, consider the following double surface integral over a pair of triangles that would appear in a Galerkin BEM 
for the Helmholtz equation: 


I = 


/xGR 


exp (ikr) 
Anr 


dS {x) dS (x), 


(34) 


where r = |x — x'|. Assume the two triangles are small and close to each other (i.e., kr < 1, see Fig.[^. The singularity 
subtraction technique works by separating the Helmholtz kernel into two pieces, a singular part (which happens to be 
the Laplace kernel) and a regular part. This is done by expanding the numerator in the Helmholtz equation as a Taylor 
series expansion around r = 0: 


I = 



1 

47rr 


dS (x') dS (x), 


(35) 


1= f f ^dSix')dS{x)+ f f \ s{x')dS{x). (36) 

The integral on the left is of the Laplace kernel, which this paper describes how to compute accurately. The integral 
on the right is completely regular, and so can be computed via standard numerical means. Because we assumed that 
kr < 1, the Taylor series is quickly convergent, and for any given accuracy, only a few terms need to be kept. 


5 Zero-Touch Case 

In the zero-touch case, because the two triangles do not touch, the double surface integral is regular and can be 
computed via standard numerical or semi-analytical means. However, in this section, we present an analytical method 
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Figure 4. The spherical coordinate system. 


for computing the integral in this case. This method uses spherical harmonics and multipole and local expansions and 
translations. Similar methods were presented in ll^l29l as part of solvers for problems in elastostatics. However, we 
build on the methods described in these references by: (1) adapting them to the kernels considered in Secs.j^an d|3}(2) 
computing the multipole expansion coefficients for a triangle exactly using Gaussian quadrature; and (3) controlling the 
error by adaptively truncating the multipole expansions and/or subdividing the problem when necessary. 


5.1 Spherical Harmonics 

The Laplace equation’s Green’s function can be expanded as 


G (x — x') = 


1 


47 r X — X' 




„ 2n + 1 r;+’ " 

n—0m——n 




(37) 


where (r', 6', (p') and (r, 0, p) are the spherical coordinates of x' and x, respectively, r< = min (r', r), and r> = 
max (r', r) (see Fig.|^. The spherical harmonics are given by OOl 

Yn = (-1)”" ^ j^j)! ) -PI™'(cos (6»)) exp , (38) 

where P™ are the associated Legendre polynomials. 

This expression can be used to build multipole and local expansions. For example, suppose we want to compute the 
potential at x due to a point source at x'. When r > r', we can build a multipole expansion: 


1 

47r |x — x'l 


^ 5] Rr{r',0\P')S::{r,e,cP). 

n—0 m——n 


(39) 


Here, 

S': (X) = s: (V, S.4.) = (», « (41) 

are the local and multipole expansion basis functions, respectively. Instead of centering the expansion around the origin, 
we can center the expansion around x*. When |x — x*| > |x' — x*|. 


1 

47r |x — x'l 


oo n 

^ ^ pr(x'-x*)5;r(x-x*). 

n=0 m— — n 


(42) 
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Figure 5. A diagram showing how to compute the double surface integral for two triangles that do not touch. 


Likewise, when |x — x*| < |x' — x*|, we can build a local expansion: 


1 

47r |x — x'l 


oo n 

n=0 m——n 


(43) 


These expressions can be used to build multipole and local expansions for arbitrary source distributions. For 
example, let us build a multipole expansion for the source distribution, p (x'), contained entirely inside an imaginary 
sphere of radius, r*, centered around x*. For a point, x, outside the imaginary sphere, the potential due to this source 
distribution is given by 


$(x) 



|<r* 


P(x') 

47r |x — x'l 


dV (x'). 


(44) 


Plugging in Eq. (|42|) and rearranging. 


E <^r(x-x*), (45) 

n=0 m——n 


where 


= 


'|x'—x*|<r* 


p(x')i?r (x'-x*)dU(xO. 


(46) 


We can build a local expansion using the same procedure. Consider a different source distribution contained entirely 
outside the imaginary sphere. For a point, x, inside the imaginary sphere, the potential due to this source distribution is 
given by 

OO n 

‘^W = E E <^«(x-x*), (47) 

n—0 m——n 


where 


'|x'—x*|>r* 


p(x') 5 r (x'-x*)dU(x'). 


(48) 


5.2 Analytical Method 

We want to compute the following double surface integral over a source triangle, S, and a receiver triangle, R: 


^ = f (cto + P • x) /" (tTg + p' • x') G (x - x') dS (x') dS (x). (49) 

Jx'^S 

First, expand the Green’s function as a multipole expansion: 


(x' - X*) (x - X*) j dS (x') dS (x), (50) 

where the expansion center, x*, is near the source triangle. Ideally, x* should be chosen so that the sphere centered 
around x* and completely containing the source triangle can be made as small as possible. There are actually only four 
possible choices of x*: the midpoints of the three edges of the source triangle and the center of the source triangle’s 
circumsphere. The one corresponding to the smallest sphere that completely contains the source triangle is chosen. 


1 = 


(cto + p • X 


/xeR 


>/ 

J x'GS 


(o-Q + p' • x') 


L n'—Q m' — — n' 
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Second, rearrange Eq. (501 by moving the double sum and the outside the inside integral; 


/ = 


f (cro+P-x)[^ ^ ff (cTp+p'-xO-R™'* (x'-x*)dS'(x')) (x-x*) j dS'(x). 

JxeR \n'=0Tn' = -n' V^x'GS / / 


(51) 

The integral over the source triangle now computes the expansion coefficients for the multipole expansion that represents 
the potential due to the linear source distribution over the source triangle: 


I = 


/xGR 


(cto + P-x) ^ ^ a;?'S';jl'(x-x*) dS'(x), 


L n'—O m' — — n' 


where 


= [ (a' + p' • x') R^:* (x' - X*) dS (xO . 

t/x'GS 


(52) 


(53) 


There are many analytical expressions available for computing the integral in Eq. ( [5^ , including those in BTII^ . In 
1 ^ . the authors presented a recursive algorithm for computing the expansion coefficients: only a™, for lower order 
and degree need to be computed explicity; the others can be computed recursively from them. However, we make 
the following observation: i?™ is polynomial, so the integrand in Eq. (53 i is polynomial. Thus, the integral can be 
computed exactly via Gaussian quadrature. A similar approach was used in ll34llT5l . although the integration domains 
in these references were lines and boxes, not triangles. We use the techniques given in for performing Gaussian 
quadrature over the triangles. 

Third, translate the multipole expansion centered around x* to a local expansion centered around y*, where y* is 
near the receiver triangle 


oo n 


I = 


/xGR 


(cro + P-x) ^ b^R^{x-y*)\dS{x). 


(54) 


V,n=0 m——n 


Ideally, y* should be chosen so that the sphere centered around y* and completely containing the receiver triangle can 
be made as small as possible. The procedure for choosing y* is the same as for choosing x*. 

Eourth, rearrange Eq. (j5^ by moving the double sum and the expansion coefficients outside the integral; 


where 


OO n / p \ 

^ = E E c / (^o+P-x)i?;r(x-y*)d5(x) , 

n—Om——n V'^xGR / 

oo n 

^ = E E 

n=0 m——n 

[ {ao + p ■ x) R^ {x - y*) dS {x). 

J xGR 


= 


(55) 


(56) 


(57) 


Like before, since ii™ is polynomial, the integral in Eq. (571 can be computed exactly via Gaussian quadrature. 

5.3 Error Control 

The expressions derived in Secs. |5 .1 1 and |5.2| involving multipole and local expansions must be truncated so that they 
can be implemented in code. Eor example, the last expression given in Sec. 5.2 becomes 


p—1 n 

^«4 = E E 

n—0 m——n 


(58) 


where only terms have been kept. Obviously, the expression is exact as p —>■ oo, but there are truncation errors 
when p < oo. These errors come from two sources: (1) the construction of the multipole expansion at the source 
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Figure 6. The radii and distances used to compute the geometric factor, p, given in Eq. (601. 


triangle; and (2) the translation of the multipole expansion to a local expansion at the receiver triangle. Luckily, these 
expansions converge rapidly in the case of the Laplace equation, so these errors can be precisely controlled by picking 
an appropriate value of p. 

Theoretical bounds as a function of p for these errors have been derived by many authors over the years. A good 
overview is given in E). The relative error is bounded by; 


e = 


Ip-I 


<Ar,P, 


where the error constant. A, depends on the problem being solved, and 

max(ri,r2) 
d- min(ri,r2)’ 


V = 


(59) 


(60) 


where ri is the radius of the multipole expansion’s bounding sphere, r2 is the radius of the local expansion’s bounding 
sphere, and d = |y* — x*| (see Fig.|^. Using Eq. (591, we can easily pick a truncation number that gives us a desired 
accuracy; 


p = 


log I ^ 1 / log (p) 


( 61 ) 


where [a;] is the ceiling of x. 

There are two issues. First, we may not always be able to use as high a value of p as we want. Since the number of 
terms in Eq. (58 < grows as p^, the memory and computational costs grow as p^ as well. Therefore, for practical reasons. 


p must be capped. In the event that p needs to be higher than this cap, we divide the larger triangle and recurse. We do 
the same thing in the event that the two bounding spheres overlap. 

Second, while values of A have been derived for special cases, such as point sources, no such values have 
been derived for the case of triangles. Instead of attempting to derive such a value analytically, we computed one 
experimentally. The experiment worked as follows. We generated 10,000 pairs of randomly placed triangles, where, 
in each pair, the two triangles did not touch. For each pair, we computed I using a semi-analytical method (using a 
high-order Gaussian quadrature for the outside integral) and Ip for p = 1, 2,..., 10. We chose a value for A so that the 
experimental data satisfied Eq. ( [59] ). We ran this experiment several times, and found A to be between 0.5 and 0.6, so 
we set A = 0.6. 

We ran a second experiment to verify that A = 0.6 works well. The experiment worked as follows. We generated 
1,000 pairs of triangles like before. For each pair, we computed I using the analytical method for several different 
choices of e: 10“^, 10“^, 10“^, and 10“®. Again, we used the value of / returned by the semi-analytical method as the 
reference value. In Fig. [7] the actual errors (along the y axis) for each pair of triangles (along the x axis) are plotted for 
each choice of e (the different colored curves). In all cases, the realized errors are below the desired error bounds. 


6 One-, Two-, and Three-Touch Cases 

In the one-, two-, and three-touch cases, because the two triangles touch, the double surface integrals are not regular, so 
they cannot be computed using standard numerical or semi-analytical means, or even the analytical method described in 
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Figure 7. This graph shows the actual relative errors when computing the double surface integrals over 1,000 randomly 
generated pairs of triangles using the analytical method for several different error bounds. In all cases, these errors were 
below the error bounds. 


Sec. 5.2 In this section, we present a method for dealing with these integrals. The approach relies on several scaling 
properties of the integrals and the kernels being integrated. The method works in roughly the following way: (1) the 
integral is broken up into several smaller integrals; (2) some of these integrals are written in terms of the original 
integral using some simple analysis; and (3) the terms are rearranged to yield an expression for the original integral that 
only requires computing regular Integrals explicitly (all other integrals are computed implicitly). 


6.1 Preliminaries: Scaling Results 

We want to compute the following double surface integral over a source triangle, S, and a receiver triangle, R: 


I = f (cto + P • x) / (cTg + p' • x') F (x - x') dS (x') dS (x), ( 62 ) 

./x^R J-x'^S 

where F (r) is any kernel that has the following scaling property: 

F (ar) = s (a) F (r). (63) 

To begin, let us break the integral into four smaller integrals by expanding the product of the two linear functions: 

I=f f (crocro + cro(p'-x') + (p-x)CTo + (p-x)(p'-x'))F(x-x')dS'(x')d 5 '(x), ( 64 ) 

J x^R J x'GS 

/ = (ToCToT^ + CTo/P +CTg/P + /PP, ( 65 ) 

where 

I^= [ f F (x - x') dS (x') dS (x), ( 66 ) 

J xGR ^ x'GS 
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Figure 8. A diagram showing the process for transforming the pair of triangles, S and R, into the pair of triangles, S* 
and R*. 


/p' = / f (p' . x') U (x - x') dS (x') dS (x), (67) 

J xGR 

/P = / f (p • x) F (x — x') dS (x') dS (x), (68) 

J xGR x^gS 

= f f (p' • x') (p • x) F (x — x') dS (x') dS (x). (69) 

J xGR 'J x^GS 

Consider two different triangles, S* and R*, which, taken together, are scaled and translated versions of S and R, 
also taken together (see Fig.[^. In other words, given a pair of points, x', x S S U R, there is a corresponding parr of 
points, y', y S S* U R*, given by 

y' = a (x' + t), y = a (x + t). (70) 

Suppose we want to compute the same integral as before, except over S* and R*; 

I*= [ [ {ao + p-y){aQ + p'■y')F{y-y')dS{y')dS{y). (71) 

The integrand is exactly the same as before. The only thing we are changing is the integration domain from S x R to 

S* X R*. Like before, let us break the integral into four smaller integrals: 

I* = + <yolf + (72) 

where 

Il= [ [ F(y-y')dF(y')dF(y), (73) 

If = [ f ip'■y')Fiy-y')dS{y')dS{y), (74) 

J yGR* *1 y'GS,K 

If=[ f {p-y)F{y-y')dS{y')dS{y), (75) 

'I y€R* ^ y'GS* 

If^=l f ip'■y'){p-y)F{y-y')dSiy')dS{y). (76) 

'J y^Rx. Jy' GS* 

Theorem 1. The four integrals over S* and R*, ll. If , If, and If can be expressed in terms of the four integrals 
over S and R, /p', /p, and /p'p: 

Il=sia)a^l\ (77) 

If = s (a) + (p' • t) I^^ , (78) 
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/P = s (a) (/P + (p • t) I^) , 

/P'P = s (a) a6 I'jp'p + (p . t) IP' + (p' • t) JP + (p • t) (p' • t) /i) . 


(79) 

(80) 


Proof. To prove Eqs. ( |77| l - ( [80] l in Theorem we are going to use the following procedure: (1) make the change of 
variables in Eq. ( [70| ; (2) use the scaling property of F (r) in Eq. (63 i; and (3) break the resulting integral into one or 
more integrals that are equal to the original four integrals, J^, /p , /p, and /p p. Eirst, let us prove Eq. (77 1 : 




F{y- y') dS (y') dS (y), 


( 81 ) 


ll = f I F {a{x + t) — a{x'+ t))dS {a{K'+ t))dS {a{x + t)), (82) 

ll = a'^ f j F {a{x — x')) dS (x') dS (x), (83) 

J x£R J x' GS 

ll = s (a) [ [ F (x — x') dS (x') dS (x), (84) 

7x6R 7x'GS 

ll = s (a) a'^. (85) 

Second, let us prove Eq. (78 1 : 

If = [ f ip'■y')Fiy-y')dS{y')dSiy), (86) 

y€R* Jy' 

If =s{a)a^ f j (p'• (x'+ t)) F (x — x') dF (x') dS” (x), (87) 

J x£R J x'GS 

/P'=s(a)a5(/P' + (pFt)/i) . (88) 

Third, let us prove Eq. (79 1 : 

If=f [ {p-y)F{y-y')dSiy')dS{y). (89) 

'J y€R* Jy' 

The analysis is exactly the same as for /p , so 

/P = s(a)a5(/P + (p.t)/i). (90) 

Eourth, let us prove Eq. (80 1 : 

If'^=l [ ip-y){p'■y')F{y-y')dSiy')dSiy), (91) 

y€R:4. ^y' GS* 

/p'p = S (a) a® / f (p • (x + t)) (p' • (x' + t)) F (x - x') dS (x') dS (x), (92) 

J x^R Jx'GS 

/p'p = s (a) (/p'p + (p • t) IP' + (p' • t) /P + (p • t) (p' • t) /i) . (93) 

At this point, we have proved Theorem 1. □ 
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Figure 9. A diagram showing how to compute the double surface integral for two triangles that share a vertex. The 
source triangle, S, is formed by the three vertices, pi, p 2 , Pa, and the receiver triangle, R, is formed by the three 
vertices, qi, qa, and qa. The two triangles share a vertex; pi = qi- 


6.2 One-Touch Case 


Consider the one-touch case in Fig.|^ Without loss of generality, assume the vertex that the two triangles share is 
located at the origin (i.e., pi = 0). The fact that they share a vertex complicates the double surface integral over them. 
We solve this problem by dividing each triangle into four smaller, congruent triangles: the source triangle, S, is divided 
into 1', 2', 3', and 4', and the receiver triangle, R, is divided into 1, 2, 3, and 4. Then, we break the integral into seven 
smaller integrals over these triangles. 

Let us look at 


/I - T^, 

1 — ± 1,1 




■ ll'Z 


ll'A 


+ ^2'R + + 


1 

3'R 


T 4 /R. 


(94) 


The subscripts denote the surfaces of integration. For example, is the integral over 1' and 1. Similarly, is the 
integral over 2' and R. The without a subscript is the original integral over S and R. In Eq. (94 1 , the six integrals, ll, 2 , 
7 ^, 3 , l\i^, / 3 /R, and correspond to pairs of triangles that do not touch. Because we know how to compute 

them, let us combine them into a single integral: 


/I 


= III 


rl 

^ remainder i 


(95) 


where 


-^remainder — -^1'2 + + ^1'4 + ^2'R + -^3'R + ^4' 


R- 


(96) 


The integral, however, has the same problem as the original integral; V and 1 share a vertex. Fortunately, V and 1 
are scaled and translated versions of S and R. We derived earlier that 


ll = s (a) . 

Since the pair of triangles, 1' and 1, is 1/2 the size of the original pair, S and R, a = 1/2, so 

'1\ 1 


A'l — ® 


2/16 




Plugging this into Eq. ( |95] l and rearranging, 

/I = f 1 - s 


— 

2/16 


-1 


n 

-‘remainder’ 


The integral, is computed implicitly in this expression. 
Now, let us look at /p ; 


(97) 


(98) 


(99) 


JP' = JP, + IP -I- 

1 -I- 1^,2 T 43,3 




-^^'R + ^3'R + ^4'R' 


( 100 ) 
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Figure 10. A diagram showing how to compute the double surface integral for two triangles that share an edge. The 
source triangle, S, is formed by the three vertices, pi, p 2 , pa, and the receiver triangle, R, is formed by the three 
vertices, qi, qa, and qa. The two triangles share an edge: p 2 = qa and pa = q 2 . 


Again, the six integrals other than are all regular, so let us combine them into a single integral: 

rp' _ rp' rp' 

1'1 ' remainder' 


( 101 ) 


This leaves We derived earlier that 


If = s (a) + (p' • t) . 


( 102 ) 


Like before, a = 1/2, but we still need to determine t. When a = 1/2, t is the point that does not change during the 
scaling and translation. This is simply the vertex that the two triangles share, so t = pi. Because we assumed that 
Pi = 0, t = 0 as well, so 


Plugging this into Eq. (|101|l and rearranging. 




(103) 


/P' = 




-1 


/P' 

remainder' 


(104) 


These same procedures can be used to compute /p and Jp P: 


/P = 1 - 


1\ 1 


2/32 


-1 


/P 

remainder’ 


/p'p = f 1 - s ( 


2/64 


remainder' 


(105) 

(106) 


6.3 Two-Touch Case 


Consider the two-touch case in Fig. Without loss of generality, assume the midpoint of the edge that the two 
triangles share is located at the origin (i.e., P 2 + Pa = 0). The fact that they share an edge complicates the double 
surface integral over them. Like in the one-touch case, we solve this problem by breaking the integral into several 
smaller integrals. 

Let us look at 




I'R + ^2'2 


■ -^213 


^ 2'4 4 ” '^3'2 4 " 4 “ 4 “ ^412 4 “ ^4/3 + ^4/4 4 “ I^i — 


(107) 
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The integrals, and are zero-touch integrals, and the integrals, l\, 2 , ll^, l\i^, 13 / 4 , /|/ 2 , -^ 4 / 3 , and l\,^, are 

one-touch integrals. Since we know how to compute them, let us combine them into a single integral: 

^ ^ 2’3 ^‘i'2 -^remainder' (108) 

That leaves J 2/3 and 13 , 2 , which have the same problem as the original integral: they correspond to parrs of triangles 
that share an edge. Fortunately, the pair, 2' and 3, and the pair, 3' and 2, are each scaled and translated versions of the 
original pair, S and R. We derived earlier that 


ll = s (a) 


(109) 


So, 


1\ 1 


ih3 = si-]-i\ ih 2 = si-]—i\ 


2/16 


3'2 


1\ 1 


Plugging these into Eq. (|108|l and rearranging. 




2/16 


n 

■^remainder' 


( 110 ) 


( 111 ) 


The integrals, 72/3 and 13 , 2 , are computed implicitly in this expression. 

Now, let us look at /p : 

-f I 212 + /^2*'3 + ^2'i + ^3'2 + ^3'3 + + -^4'2 + ^4'3 + + ^^Sl ~ 

Again, let us combine the ten integrals that we know how to compute into one: 


fP = fP -4- fP + fP 
^ ^2'3 ' ^3'2 ^ '^remainder’ 


We derived earlier that 


So, 


if = s (a) (^/P' -f (p' • t) . 

I2'3 + l3'2 = S (^) (0 ^ (P' ■ + t3'2)) /'■ 


( 112 ) 

(113) 

(114) 

(115) 


We need to determine t 2'3 and ty 2 . As discussed earlier, when a = 1/2, the translation vector, t, corresponds to the 
point that does not change during the transformation. Thus, for the pair, 2' and 3, that is p 2 , and for the pair, 2 and 3', 
that is P 3 . The edge that S and R share is centered around the origin, so t 2'3 -f 13/2 = 0, which means that 


/p„ + /p;- s( 


2/16 


Plugging this into Eq. (|1 13 1 and rearranging. 


1\ 1 


/P = 1 - s U /.P 


The same analysis can be used for comuting /P: 


/P = 1 - s 


2/16 


1\ 1 


2/16 


^ remainder' 


remainder' 


(116) 


(117) 


(118) 


Einally, let us look at /p P: 


where 


(119) 

rP'p _ tp'p I 7P'P 1 rP'p I rP'p I 7P'P I rP'p _i_ rp'p , jp'p , jp'p _ rp'p 

-'remainder ~ -'I'R ^2'2 ~^^2'4 ^3'3 ~^-^3'i ^4'2 "I"-'4'3 ~^^4’4 ^S1 -'I'l ’ 
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Figure 11. A diagram showing how to compute the double surface integral for two triangles that are the same. The 
source triangle, S, is formed by the three vertices, pi, p 2 , pa, and the receiver triangle, R, is formed by the three 
vertices, qi, q 2 , and qa. The two triangles are the same: pi = qi, p 2 = qa, and pa = qa- 


We derived earlier that 


/P'P = s (a) a® (/P'P + (p • t) /P' + (p' • t) JP + (p • t) (p' • t) /i) . 


( 121 ) 


So, 


= s 


The terms linear in p' and p disappear because t 2 'a + t 3'2 = 0, but unfortunately, the other terms remain. Plugging Eq. 
\\22\ into Eq. ( |1 19| l and rearranging, 

-1 


/PP= (l-s( M 


2/32 


( rp'p „p'p^ 

l remainder ' J ’ 


where 


1\ 1 


a*" = s ( ^ j ^ ((P • t 2 'a) (p' ■ t 2 'a) + (p • ta' 2 ) (p' • ta' 2 ))/^ 


(123) 


(124) 


6.4 Three-Touch Case 

Consider the three-touch case in Eig. [TT] In this case, the two triangles are the same. Without loss of generality, assume 
that the centroids of the two triangles are located at the origin (i.e., Pi + P 2 + Pa = 0)- The fact that the two triangles 
are the same complicates the double surface integral over them. Like in the one- and two-touch cases, we solve this 
problem by breaking the integral into several smaller integrals. 

Let us look at 


/I = I}r 


The integrals, ll, 2 , Ii/a, ^a'l’ T 3 / 2 . are one-touch integrals, and the integrals, ll,^, I 214 , 1 ^/ 4 , /j/a, /|/ 2 ’ 

and l\,^, are two-touch integrals. Like before, because we know how to compute them, let us combine these 12 integrals 
into a single integral: 

= ll'l + 1^2 + Ta'a + (126) 

That leaves the four integrals, 1 ^, 2 , l\i^, and l\n, which have the same problem as the original: they correspond to 
pairs of triangles that are the same. We derived earlier that 

Il=s{a)a^l\ (127) 


So, 


Il'i 


+ 4 


2'2 


43 = 
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The pair of triangles, 4' and 4, are not only scaled by a factor of 1/2, but also rotated by 180 degrees. Because they are 
centered around the origin, we can achieve this rotation by setting a = —1/2: 


L 


1 

4'4 


= S 




Plugging these into Eq. \\26\ and rearranging. 


/! = 





ri 

remainder ■ 


(129) 


(130) 


The same analysis from here and before can be used to compute expressions for /p , Jp, and /p P: 

-1 


/P = 1 - 


1\ 3 


2/32 


1\ 3 


+ s - 


1 


/P= 1-s IL 


2/32 


2/32 


1\ 1 


/P 

remainder’ 


1\ 3 


/PP= 1 -s X -xhrr 


2/64 


1\ 1 


2/32 
-1 


-1 


^remainder’ 


2/64 


(itL 


remainder 


tP P 


where 




and ti>i = pi, t2'2 = P2, and tya = pa. 


(131) 

(132) 

(133) 

(134) 


7 Numerical Examples 

We have implemented the two methods that were described in Secs.|^and|^in MATLAB as part of a Galerkin BEM 
library for the Laplace equation. Currently, the library only supports external problems with Dirichlet boundary 
conditions, but even such simple problems provide an opportunity for demonstrating these two methods. The library 
supports constant and linear triangular elements, and for comparison purposes, we have also implemented the collocation 
method. Thus, the library contains four solvers: (1) constant collocation; (2) linear collocation; (3) constant Galerkin; 
and (4) linear Galerkin. In this section, we demonstrate that the two methods work as advertised by using the library 
to solve two test problems. To begin, we mention a couple implementation details of the library. Eollowing that, we 
describe and give the analytical solutions to the test problems, and then show some results from solving the problems 
using the library. 


7.1 Implementation Details 

There are two implementation details that should be mentioned. 

Eirst, in MATLAB, for the same accuracy, the analytical method for the zero-touch case described in Sec.|^is faster 
for parrs of triangles that are far away from each other, and the semi-analytical method is faster for pairs of triangles 
that are nearby each other. Thus, when a zero-touch integral needs to be computed directly, and not indirectly as part of 
a one-, two-, or three-touch integral, the analytical method is used. However, the semi-analytical method is used instead 
when a zero-touch integral needs to be computed during the computation of a one-, two-, or three-touch integral. Eor 
example, when computing a one-touch integral, six zero-touch integrals need to be computed, and these are done so 
semi-analytically. 

Second, for external problems with Dirichlet boundary conditions, the only kernel that needs to be integrated is the 
Laplace equation’s Green’s function. That is, when F (r) = G (r), where 


G(r) 


1 

47 r |r| 


(135) 
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Potential, Analytical Solution Potential, Analytical Solution 



X X 

Figure 12. The potential due to a point source placed inside a cube (left) or a sphere (right). 


The associated scaling function, s (a), for this kernel is given by 

1 


I (a) = 


a 


Below, we quickly specialize and summarize the expressions derived in Sec.j^for this kernel. 
For the one-touch case; 


Q 

ft _ ° rl jp' — mjP jp _ djP jp'p _ ^jP P 


= ¥■ 


remainder ’ 


16 

remainder ’ 


16 

“^remainder’ 


32 

“^remainder’ 


For the two-touch case: 


8 


16 


rl _ Irl TP' _ ^rP tP — ^ tP tP P _ 

remainder’ y remainder ’ 


{i± 


remainder 


where 




For the three-touch case: 


ri = ori TP' = -J'P' f P = - fP 

remainder’ ^ rr'^remainder’ ^7 rei 


^ remainder ’ 


JP P _ _ 


{ti. 


remainder 


- 1 - 0 ^ 


where 


aPP = 


1 


((P ■ tl'l) (P' ■ tl'l) -P (p • t2'2) (P' • t2'2) + (P • ts's) (P' ■ ts's)) 


(136) 

(137) 

(138) 

(139) 

(140) 

(141) 


7.2 Test Problem #1: Cube Behaving Like a Point Source 


Consider a cube, C, centered around the origin, and place a point source at Xq, where Xq is inside the cube. The 
potential due to the point source is given by 


0’’' (x) 


1 

47r |x - Xo| 


(142) 


and is shown in Fig. Suppose we want to compute an equivalent source distribution on the surface of the cube that 
gives rise to the same potential outside the cube as that of the pont source. In other words, we want to compute a surface 
source distribution on the cube such that 


cj) {x) = {x) ^ -r, x&dC, 

Itt |x - xqI 

where dC is the surface of the cube. Obviously, the solution is (p (x) = cjf^ (x). 


(143) 
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Figure 13. Two meshes used to compute the numerical solutions to the two test problems. On the left: a mesh for the 
cube (300 triangles). On the right: a mesh for the sphere (204 triangles). 


7.3 Test Problem #2: Sphere Behaving Like a Point Source 

The sphere test problem is the same as the cube test problem, except that the boundary is a sphere instead of a cube. In 
other words, we want to compute a source distribution on the surface of the sphere that gives rise to a potential outside 
of the sphere that is the same as that due to a point source placed inside the sphere (see Fig.[T2|). The point source’s 
position is the same in both test problems. The only difference is the geometry of the boundary. 


7.4 Error Analysis 


First, we had to construct triangular meshes for the two test problems. Fig.[T^shows two example meshes, one for each 
test problem. The meshes were constructed so that all the triangles were as close to equilateral and the same size as 
each other as possible. In fact, in the cube mesh, all the triangles are congruent. 

Second, we solved the two test problems using the four solvers. The software performed very well for the cube 


test problem. In particular, the linear Galerkin solver performed extremely well. Fig. 14 shows the relative error in the 
volume around the cube for each solver when the mesh had 300 triangles. To show that the relative errors decrease 
as the mesh size increases, we computed and plotted the maximum relative error in the volume around the cube as a 
function of the number of unknowns, N (see Fig.[T6)l. For constant triangular elements, the number of unknowns is the 
number of triangles. For linear triangular elements, the number of unknowns is the number of vertices. All four curves 
behave at least as good as 1/N. In general, linear triangular elements are better than constant triangular elements, and 
the Galerkin methods are better than the collocation methods. The linear Galerkin method behaves even better than 
1 /N for larger N. The nice thing about the cube test problem is that the boundary is piecewise planar, and so can be 
modeled exactly using planar triangular elements. Thus, there are no errors from approximating the geometry of the 
boundary. 

The software performed very well for the sphere test problem as well. Fig. [TS] shows the relative error in the volume 
around the sphere for each solver when the mesh had 204 triangles. Fig. |16| shows the maximum relative error as a 
function of the number of unknowns. Like the cube test problem, all four curves behave at least as good asas 1/A^, and 
the linear Galerkin method is the best, decaying nearly as 1 /N"^ for larger N. 


8 Conclusion 

We have presented a method for computing the double surface integrals encountered in the Galerkin BEM. When the 
boundary is discretized using triangular elements, these integrals are performed over pairs of these triangles. They can 
be extremely difficult to compute, especially when the two triangles share a vertex, an edge, or are the same. This is 
because the kernels being integrated are often singular along the corners and edges of these triangles. We have solved 
this problem by using several scaling properties of the integrals and the kernels being integrated. The integral is broken 
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Figure 14. The relative error of the numerical solutions returned by the four solvers for the cube test problem using the 
mesh seen in Fig. 13 The value for N shown in the titles is the number of unknowns. For constant elements, that is the 
number of triangles. For linear elements, that is the number of vertices. 


up into several smaller ones, some of which are written in terms of the original. This is done in such a way that only 
completely regular integrals have to be computed explicitly. 

We have also presented an analytical method for computing the integrals when the two triangles do not touch. The 
method uses spherical harmonics and multipole and local expansions and translations. The only source of error in this 
method is how soon to truncate these expansions. However, the truncation number is adaptively selected to achieve a 
desired error bound. 

Finally, we have implemented the two methods as part of a Galerkin BEM library using MATLAB, and we have 
made this library freely available for download from our webpage. 
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Figure 15. The relative error of the numerical solutions returned by the four solvers for the sphere test problem using 


the mesh seen in Fig. 13 
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